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ABSTRACT 

We present simulations of the non-linear evolution of streaming instabilities in protoplanetary disks. 
The two components of the disk, gas treated with grid hydrodynamics and solids treated as super- 
particles, are mutually coupled by drag forces. We find that the initially laminar equilibrium flow 
spontaneously develops into turbulence in our unstratified local model. Marginally coupled solids 
(that couple to the gas on a Keplerian time-scale) trigger an upward cascade to large particle clumps 
with peak overdensities above 100. The clumps evolve dynamically by losing material downstream to 
the radial drift flow while receiving recycled material from upstream. Smaller, more tightly coupled 
solids produce weaker turbulence with more transient overdensities on smaller length scales. The net 
inward radial drift is decreased for marginally coupled particles, whereas the tightly coupled particles 
migrate faster in the saturated turbulent state. The turbulent diffusion of solid particles, measured 
by their random walk, depends strongly on their stopping time and on the solids-to-gas ratio of the 
background state, but diffusion is generally modest, particularly for tightly coupled solids. Angular 
momentum transport is too weak and of the wrong sign to influence stellar accretion. Self-gravity and 
collisions will be needed to determine the relevance of particle overdensities for planetesimal formation. 
Subject headings: diffusion — hydrodynamics — instabilities — planetary systems: protoplanetary 
disks — solar system: formation — turbulence 



L INTRODUCTION 

This paper extends our preparato ry numerical investi- 
gatio ns of the streaming instability (jYoudin fc JohansenI 
120071 hereafter YJ) into the non-linear regime. The 
linear streaming in s tabili ty was first described by 
lYoudin fc GoodmanI ([2005', hereafter YG) who found 
that the radial and azimuthal drift of solids through 
gas in a protoplanetary disk triggers growing oscillations 
that concentrate particles. YJ details the numerical tech- 
niques used to study the evolution of solids and gas in 
a local patch of a protoplanetary disk and demonstrates 
that our code successfully reproduces the linear growth 
rates derived by YG. This paper describes the non-linear 
evolution of the streaming instability to a fully turbulent 
state and studies the consequences for particle concen- 
tration and transport. 

The starward drift of solids, caused by the sub- 
Keplerian headwind encountered by the particles, is not 
just a trigger for streaming instabilities, but also a source 
of theoretical difficulties. Growing planetesimals by co- 
agulation faces severe time-scale constraints due to the 
loss off solids (ultimately to the star or sublimation in 
the inner disk). The restriction is most acute for 10 cm 
"rocks" through 1 m "boulders" with drift times of only a 
few hundred orbits (|Weidenschillinel[l977a[ ) in most of the 
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planet-f orming region of standard minimum m ass nebula 
models (IWeidenschilling|ll977bl : IHavashilHOSl . The drift 
of mm-sized solids in a few times 10^ years at 30 AU is at 
best marginally consistent with the observed mm-excess 
from the outer parts T-Tauri disks with ages of a few 
Mvr (IWilner et al.l [2000HRodmann et al.l 120061: see also 
iBrauer et al.l 1200^ for possible ways to maintain such a 
population of "pebbles"). This mismatch between the- 
ory and observations may indicate that simple drift time 
estimates are missing important dynamical effects. 

Several mechanisms could impede the radial influx of 
solids. The increased inertia of solids in a dense midplane 
sublayer decreases drift speeds as the local gas mass frac- 
tion squared (Nakagawa et al. 1986; Youdin & Chiang 
[200l . Since such high densi ties may trigger rapid grav- 
itational collapse of sohds (jYoudin fc Shul l2002f l. sed- 
imentation alone is not a satisfactory explanation of 
the long lifetimes of pebbles in the disks, even if the 
turbulence is weak enough to allow the formation of 
an extremely thin mid-plane layer. Turbulent diffu- 
sion in accretion disks will maintain a small fraction of 
particles in the oute r disk (jStepinski fc Valagead 119961 : 
iTakeuchi fc Lh]|2002l ). but this scenario requires a parti- 
cle reservoir that exceeds by far the observed amount 
of mm-sized solids and thus implies disk masses that 
are orders of magnitude larger th an minimum mass 
mode l s. Giant anticyclonic vortices (iBarge fc Sommerial 
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119951: Ide la Fuente Marcos fc Bared 12001 " staU migra- 
tion by trapping marginally coupled solids. How- 
ever the f ormation and stability of vortices in disks is 
not c l ear (iGoodman et all 19871: iKlahr fc Bodenheimeil 
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iFromang fc Nelsoiill2005D . Any local press ure maximum 
- no t only vortices - can trap b oulders (jKlahr fc LinI 
1200 It iHaghighipour fc Bosd l2 003l). e.g. spiral arms of 
massive self-gravitating disks ()Rice et al.ll2004[ ) or even 
transient p ressure enhancements in magnetorotational 
turbulence (| Johansen et al.l l2006b| ). The present work 
will show that streaming turbulence modestly slows the 
average radial drift of marginally coupled solids. An ul- 
timate solution of the drift problem may require rapid 
(faster than drift) planetesimal formation (by gravi- 
tational collapse and/or coagulation) and fragmenta- 
tion to maintain observed p opulations of small solids 
(jPullemond fc Dominikll2005l ). 

The dynamical particle trapping mechanisms men- 
tioned above increase particle densities, with an efficiency 
that depends on (often uncertain) structure lifetimes. 
Local particle overdensities can seed gravitational col- 
lapse of solids and affect the rates of (and balance be- 
tween) coagulation and coUisional fragmentation. Opti- 
cally thick clumps could even influence radiative transfer 
if the disk itself is optically thin, and thereby alter ob- 
servational e stimates of disk mass and particle size (see 
iDraind [20061 for a general discussion of the role of op- 
tical depth, but not clumping per se). Radial drift in- 
herently augments the surface density of solids in the 
inner disk as particles pile up from larger orbital radii, 
as long as particles are smaller tha ,n the gas mean free 
path so that Epstein d rag applies (|Youdin fc Shul[200l : 
I Youdin fc Chian^l2004f l. In simulations and experiments 



of forced Kolmogorov turbulence, particles concentrate 
in low vorticity reg i ons at the viscous d issipation scale 
(|Fessler et al.H'l994|- iPadoan et al.ll2006f) . Efficient col- 
lection requires small particles that couple to the rapid 
turnover time at the dissipation scale. ICuzzi et all ()2001l ) 
apply this passive turbulent concentration to the size- 
sorting of chondrules (abundant, partially-molten, mm- 
sized inclusions found in primitive meteorites) in the in- 
n er solar nebula. 

iJohansen et al.l (|2006aD discovered active turbulent 
concentration (active meaning that the drag feedback 
on gas was included) of larger particles (from cm-sized 
pebbles to m-sized boulders) in simulations of Kclvin- 
Helmholz midplane turbulence. Dense clumps of solids 
plow through the gas at near the Keplerian speed, scoop- 
ing up more isolated particles that move with the sub- 
Keplerian headwind. Since this particle concentration 
mechanism relies on two-way drag forces, it was (and still 
is) considered a non-linear manifestation of streaming in- 
stabilities. The current work further explores active con- 
centration, isolating the role of drag feedback by ignoring 
vertical stratification. We consider different geometries, 
both axisy mmetric in the r adial- v ertical plane and fully 
3-D, than I Johansen et"aLl (|2006a| ). who considered the 
azimuthal- vertical plane, and we use the higher order in- 
terpolation scheme for drag forces described in YJ. We 
will thus show that the "pure" streaming instability also 
produces strongly non-linear particle overdensities. 

Turbulent diffusion controls the extent to which par- 
ticles sediment in the mid-plane (|Dubrulle et al.|[T995D 
and whether (or how fast) self- gravity can coll ect solids 
into rings and bound clumps (jYoudinI l2005al |H). Dif- 
fusion is the most fundamental parameter governing 
whether planetesimals can form by gravitational in- 
stability (as originally proposed by ISafronovl 1X9691 : 
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Fig. 1. — Peak growth rate s of the streaming instability versus 
the solids-to-gas ratio e for friction times of = 1.0 (solid line) 
and Tb = 0.1 (dashed line). The steep rise in growth rate when 
Ts = 0.1 particles cross Pp/pg = 1 explains the cavitation in the 
non-linear run AB (see Fig. [5}. 



iGoldreich fc Wardll973^ ■ because the oft-mentioned crit- 
ical density for gravitational collapse is irrelevant when 
drag forces are included to transf er angular m omentum 
from the solids to the gas (Wardl ll976[ l2000f ). Passive 
diffusion of particles in magnetorotatio nal turbulence 
has been found to be quite strong (.Johansen fc Klahil 
2005t lTurne_ret all [20061 : iFromang fc Papaloizoul l2006t 



CarbalHdo et al.l [20061 ). although the Schmidt number 



(the ratio between the turbulent viscosity of the gas 
and the particle diffusion) increases i n the presence 
of a net vertical mag netic field (Carballido et al.l 120051 : 
I Johansen et al.l[2b06cl ). In the present work we measure 
the active particle diffusion in streaming turbulence, by 
considering the random walk of the particles away from 
a reference point, and find it to be relatively weak, espe- 
cially for smaller particles. 

The paper is built up as follows. In fj2|we briefly reit- 
erate the physical model of the protoplanetary disk and 
our numerical method for solving the dynamical equa- 
tions of gas and solids. In fj3| we present the non-linear 
simulations and the topography of the turbulent state, 
before analyzing in ^j4|the statistics and causes of particle 
clumping in more detail. Then ^ addresses the trans- 
port and diffusion of particles and angular momentum in 
the saturated streaming turbulence. We summarize our 
results in 

2. NUMERICAL METHOD 

The dynamical equations and the numerical method 
are presented in detail in YJ, but we bricffy recapitulate 
the main points in this section. As a numerical solver we 
use the Pencil Code'^ . This is a finite difference code that 
uses 6th order symmetric spatial derivati ves and a 3rd 
order Runge-Kutta time integration (see iBrandenburel 
[200l for details). 

We model a local patch in a protoplanetary 
disk with the shearing s heet approximation (e.g. 
IGoldreich fc Tremaind [19781 ). A Cartesian coordinate 
frame that corotates with the Kepler frequency J7 at a 
distance r from the central star is oriented with x, y 
and z axes pointing radially outwards, along the orbital 
direction, and vertically out of the disk (parallel to the 

^ See |http : //www. nordita. dk/software/pencll- code/ [ 
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TABLE 1 

Run Parameters 



Run 




e 




X Ny 


X 




At 


AA 


0.1 


0.2 


4.0 X 4.0 X 4.0 


256 X 


1 X 256 




2000.0 


AB 


0.1 


1.0 


2.0 X 2.0 X 2.0 


256 X 


1 X 256 


1.6 X 10"^ 


50.0 


AC 


0.1 


3.0 


2.0 X 2.0 X 2.0 


256 X 


1 X 256 


1.6 X 10® 


50.0 


BA 


1.0 


0.2 


40.0 X 40.0 X 40.0 


256 X 


1 X 256 


1.6 X 10® 


500.0 


BB 


1.0 


1.0 


20.0 X 20.0 X 20.0 


256 X 


1 X 256 


1.6 X 10® 


250.0 


BC 


1.0 


3.0 


20.0 X 20.0 X 20.0 


256 X 


1 X 256 


1.6 X 10® 


250.0 


AB-3D 


0.1 


1.0 


2.0 X 2.0 X 2.0 


128 X 128 X 128 


2.0 X lO'^ 


35.0 


BA-3D 


1.0 


0.2 


40.0 X 40.0 X 40.0 


128 X 128 X 128 


2.0 X lO'^ 


300.0 



Note. — Col. (1): Name of run. Col. (2): Friction time. Col. (3): Solids-to- 
gas ratio. Col. (4): Box size in units of r)r. Col. (5): Grid resolution. Col. (6); 
Number of particles. Col. (7): Total run time in units of f!~^ . 



Keplerian angular momentum vector), respectively. We 
solve the equations of motion for deviations from Keple- 
rian rotation in an unstratified model, i.e. vertical gravity 
is ignored. The gas (and not the sohds) is subject to pres- 
sure forces, including the global radial pressure gradient, 
constant in the local approximation and measured by the 
dimensionless parameter 



dP/dr 



Cs 



(1) 



where vk — is the Keplerian orbital speed, while P, 
Pg, and Cs are the pressure, density and sound speed of 
our isothermal gas. All our simulations use rj — 0.005 
and Cs/wK = H/r = 0.1, where H is the gas scale- height. 
Our results can be applied to different values of rj if veloc- 
ities are scaled by yywK, the pressure-supported velocity, 
and lengths are scales by ryr, the radial distance between 
points where Keplerian and pressure-supported velocities 
are equal.* 

The solids are treated alternatively as a pressureless 
fluid or as superparticles that each contain the mass of 
many actual solid bodies. Solids and gas mutually inter- 
act by frictional drag forces that are linear in the rela- 
tive velocity. This models small particles with a friction 
(or stopping) time Tf that is independent of the relative 
speed between gas and particles (i.e. no turbulent wakes 
form). The translation from friction time to the radius 
of a particle depends only on gas properties and the ma- 
terial density of the solids. As a rule of thumb the radius 
of a compact icy particle in meters is roughly equal to 
the dimensionless stopping time 



(2) 



at Jupiter's location (r » 5 AU) in standard minimum 
mass nebula models (|Havashilll981[ ). 

When solids are treated as numerical particles, we cal- 
culate the drag acceleration by interpolating the gas ve- 
locity at the positions of the particles using a second- 
order spline fit to the 9 (27) nearest grid points that 
surround a given particle for 2-D (3-D) grids. To con- 
serve momentum we assign the drag force on each sin- 
gle particle back to the ga s using a Triangular Shaped 
Cloud (TSC) scheme (Hocknev fc East wood 1981). This 
smoothing of the particles' influence helps overcome shot 

The value of H/ {rjr) (= 20 in our simulations) changes with this 
scaling, but should not affect the results as long as Mach numbers 
remain low (H/r <g 1) so that gas compressibility is insignificant. 



noise, and should not be seen as an SPH (smoothed par- 
ticle hydrodynamics) approach since our particles carry 
no hydrodynamic properties. We showed in YJ that 1 
particle per grid point is enough to resolve the linear 
growth of the streaming instability, but to better handle 
Poisson fluctuations for a wider range of densities, we 
generally use 25 particles per grid point in the non-linear 
simulations. 

An equilibrium solution to the coupled equations 
of motion of t h e gas and the solids was found by 
iNakagawa et al.l (|1986l hereafter referred to as NSH) 
where drag balances the radial pressure gradient and 
Coriolis forces. YG found that this equilibrium triggered 
a linear drag instability. The peak growth rate of this 
streaming instability is shown in Fig. [T] as a function of 
the solids-to-gas ratio e for Tg of 0.1 and 1.0. See Figs. 1 
and 2 of YJ (and the accompanying text) for the depen- 
dence of growth rates on wavenumber. 

3. NON-LINEAR EVOLUTION TO TURBULENCE 

With confidence from YJ that the code solves correctly 
for the linear growth of the streaming instability, we turn 
our focus to the non-linear evolution into turbulence. We 
generically refer to the non-linear states of our runs as 
"turbulent," because they contain stochastic fluctuations 
that diffuse material and momentum. Some cases (the 
gas-dominated A A and BA runs, see below) appear more 
wave-like with peaks in the spatial and temporal Fourier 
spectra (as we will see in Fig. [T2|) . However, even these 
runs exhibit diffusion and stochastic fluctuations on a 
range of scales, so we also label them as turbulent. This 
section describes the simulation parameters and main re- 
sults for marginal vs. tighter coupling. More detailed 
analyses of the turbulent state follow in later sections. 

3.1. Run Parameters and Initialization 

The parameters of the different non-linear simulations 
are listed in Table [T] We consider two particle sizes, 
represented as friction times: tightly coupled solids with 
Ts = 0.1 (those runs are labeled A*, where * represents 
a solids-to-gas ratio label) and larger, more loosely cou- 
pled particles with Ts = 1.0 (labeled B*). Three values 
of the sohds-to-gas mass ratio, e — 0.2,1.0,3.0 (labeled 
*A, *B, *C, respectively) are considered for each friction 
time. For instance model AB uses Tg = 0.1 and e = 1.0. 
The chosen particle abundances are well above the Solar 
composition of e ~ 0.01, but can very well be achieved 
in a sedimented mid-plane layer of solids, depending on 
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Fig. 2. — Particle density snapshots for run BA with friction time Ts = 1.0 and a solids-to-gas ratio of e = 0.2. Particle densities increase 
from black (zero density) to bright yellow/ white (solids-to-gas of unity or higher). The evident linear wavelength in the first frame results 
from the streaming instability feeding off the drift of the particles through the gas. Subsequent frames document a surprising consequence 
of the self-consistently generated turbulence: the non-linear cascade of dense particle clumps into larger filaments. 



turbulent diffusion and various particle enrichment or gas 
depletion mechanisms. 

The size of the simulation box was chosen in all cases 
such that the most unstable radial wavelengths are re- 
solved with at least 8 grid points. Two of the runs (la- 
beled AB-3D, BA-3D) were fully 3-D, aU others were 2.5- 
D simulations of the radial-vertical plane with all three 
velocity components, consistent with the linear analy- 
sis of YG. Fully periodic boundary conditions were used 
for the 2.5-D runs, while the 3-D simulations impose a 
shear-periodic boundary condition in the radial direction 
(see §3.2.1 of YJ). Particles are initially placed randomly 
throughout the simulation box. This "warm start" gives 
a white noise power spectrum with scale-independent 



Fourier amplitudes of pp{k)/{pp) 



in the par- 



ticle density. The noise serves as a seed for streaming 



instabilities. The velocities of gas and solids are initially 
set to the equilibrium values of NSH. 

3.2. Marginally Coupled Boulders 

Many drag force phenomena are most prominent for 
marginally coupled, Ts = 1, particles, corresponding to 
approximately meter-sized boulders at r « 5 AU in the 
solar nebula. Streaming instabilities are no exception, 
with fast linear growth^ and significant particle clumping 
in this regime. Fig. [2] shows four snapshots of the evolu- 
tion of the streaming instability into turbulence for run 
BA (ts = 1.0 e = 0.2). The initial growth is dominated 
by the fastest linear modes (first frame of Fig. [2]), consis- 
tent with the maximum analytic growth rate, s ~ 0.1f2 

^ Somewhat paradoxically, tight coupling gives faster growth in 
the particle-dominated regime, but on smaller scales. 
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Fig. 3. — The saturated state of runs BB and BC (both at a time of t = 100Q~^). The range of the soHds-to-gas ratio is from 0-5 in 
the left plot and from 0-15 in the right plot, giving both the same relative scale for particle overdensities as Fig. [2] The tendency for dense 
clumps to lean against the radial drift flow is evident. 



for k^rir k, 1 (see Fig. [T] and also Figs. 1 & 2 of YJ). 

A non-linearly fluctuating, i.e. turbulent, state is 
reached after some 80 local shear times (second frame 
of Fig. [2]). Solids become concentrated in a few massive 
clumps surrounded by an ocean of lower density mate- 
rial. Radial drift speeds are lower in such dense regions 
(we discuss the reduced radial drift further in §5.ip . Solid 
particles are eventually lost downstream from the clumps 
into the voids, where the radial drift is faster, until they 
fall into another dense particle clump. Over a time-scale 
of more than 100 shear times (third and fourth frame 
of Fig. (2) this leads to an upward cascade of the den- 
sity structure into extended filaments (actually rolls and 
sheets if we extend into the symmetric azimuthal dimen- 
sion). The filaments are predominantly aligned in the 
vertical direction, which maximizes their ability to in- 
tercept particles, but are slightly tilted radially in al- 
ternating directions. Strong bulk motions are exhibited 
by the filaments along their long axis. This helps them 
stay upwind (motions are in the +z, -\-x or — z, -I-.t di- 
rections), and leads to their disruption in several orbital 
times when alternately aligned filaments collide.^ The 
bulk motion also leads to efficient mixing of particles, es- 
pecially in the vertical direction (see ijS.Sp . The extended 
filaments are closely related to the long-lived vertically- 
oscillating clumps seen in 2-D simulations of the Kelvin- 
H elmholtz instabihty wi th Ts = 1.0 particles (see Fig. 8 
of lJohansen et al]l2006a| ). 

The Ts = 1.0 runs with larger e values (BB and BC) 
evolve similarly to run BA, but with a less pronounced 
cascade to larger scales (the saturated states of those 
two runs are shown in Fig. [3]) . Fig. d] shows the evo- 
lution of maximum particle density (assigned with the 
TSC scheme to the grid) versus time for all three runs. 
The non-linear state is characterized by density peaks 

® The behavior described is best seen in a movie of run AB wrhich 

can be downloaded from 

|http : //www . mpia . de/homes/ j ohansen/research.en ■ php| 
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Fig. 4. — Maximum bulk density of solids, in units of the average 
gas density in the box, as a function of time for the three marginally 
coupled runs. The maximum density is generally around two or- 
ders of magnitude higher than the mean bulk density of the solids. 
The particle density has been assigned to the mesh using the TSC 
scheme. 

100 times (or more) above the average particle density. 
Run BA has a longer run time, not only because the gas- 
dominated case is more astrophysically interesting, but 
also because it took longer to reach a saturated state. 
Fig. [4] shows signs of secular growth in densities over 
the entire At = 500J7~^. Even longer runs would bet- 
ter characterize long term fluctuations in the saturated 
state, but such computational resources would probably 
be better spent on a more realistic model with vertical 
gravity. 

Table [5] lists turbulent Mach numbers of the gas 
flow (after subtracting the mean flow, see ii4.ip . The 
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TABLE 2 
Flow Properties 



Run 


Ts 


e 






MHy 


Ma^ 




Vx 1i 


;^(NSH) 


AA 


0.1 


0.2 


5.7 X 10- 


■4 


1.1 X 10-3 


3.4 X 10" 


-3 


-0.138 


-0.138 


AB 


0.1 


1.0 


1.2 X 10" 


■2 


6.1 X 10-3 


8.5 X 10" 


-3 


-0.108 


-0.050 


AC 


0.1 


3.0 


8.7 X 10- 


3 


4.5 X 10-3 


6.4 X 10" 


-3 


-0.035 


-0.012 


BA 


1.0 


0.2 


1.2 X 10- 


-2 


1.8 X 10-2 


4.0 X 10" 


-2 


-0.520 


-0.820 


BB 


1.0 


1.0 


9.3 X 10- 


3 


1.1 X 10-2 


9.2 X 10" 


-3 


-0.341 


-0.400 


BC 


1.0 


3.0 


8.9 X 10- 


3 


1.3 X 10-2 


1.1 X 10" 


-2 


-0.118 


-0.118 


AB-3D 


0.1 


1.0 


5.3 X 10- 


3 


3.4 X 10-3 


2.7 X 10" 


-3 


-0.064 


-0.050 


BA-3D 


1.0 


0.2 


1.2 X 10- 


■2 


1.7 X 10-2 


3.3 X 10" 


-2 


-0.545 


-0.820 


Note. 


— Col. (1): 


Name of run. Col. (2): 


Friction time. 


Col. (3) 


: Solids- 



to-gas ratio. Col. (4)-(6): Turbulent Mach number of the gas. Col. (7): Mean 
radial particle velocity in units of rjv-K- Col. (8): Moan radial particle velocity in 
NSH state. 
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Fig. 5. — The onset of streaming turbulence for run AB (ts = 0.1, e = 1.0). The plots show color coded particle density (black is zero 
particle density, bright is a solids-to-gas ratio of 5 or higher). The first frame shows only the initial Poisson noise. After around one orbital 
period small voids form. The inner edges of the cavities are loaded with particles that can fall rapidly through the voids in the absence of 
any collective drag force effects there. The voids rapidly expand, and a self-sustained turbulent state sets in after around 2 orbits. This 
atypical onset of turbulence is caused by an increased growth rate of the streaming instability in slightly overdense grid cells (see Fig. [T]l. 
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anisotropic turbulence of run BA, with stronger fluctu- 
ations in the vertical direction, is clear. Turbulence is 
more isotropic in the other marginally coupled runs. 

3.3. Tightly Coupled Rocks 

Simulations with shorter friction times are more costly 
because the shorter unstable length-scales (Figs. 1 & 2 
of YJ) impose more stringent Courant criteria. The ef- 
fort was nevertheless rewarded with a qualitatively very 
different behavior for the = 0.1 runs. These particles 
correspond to solid rocks of approximately 10 cm size at 
r — 5 AU in the solar nebula. 

3.3.1. Ts = 0.1, e — 1.0; Cavitation 

Fig. [5] shows four snapshots of the particle density 
for run AB (friction time Ts — 0.1, solids-to-gas ratio 
e — 1.0). The first frame displays the initial Poisson 
noise. In contrast to run BA (see the first frame of 
Fig. [2]), we do not see the smooth growth of linear waves 
over ten or more of orbital times (an expectation which 
follows from the peak growth rate s « 0.15i7). Instead a 
few voids with dense inner rims appear hy t — 6f2~^ (sec- 
ond frame). The cavities expand rapidly (third frame), 
leading to a fully turbulent state after only two orbits, 
i.e. t ~ \2fl~^ (fourth frame). 

The effect of Poisson fluctuations on the linear growth 
properties in Fig.[T]largely explains the surprisingly rapid 
and non-uniform growth. The e = 1.0 state lies amid 
a steep rise in growth rates from the gas-dominated to 
particle-dominated regimes. Speciflcally the growth time 
for e = 1.0, tgrow = 1/s = 6.8I?~^, is halved for a modest 
increase in the solids-to-gas ratio to e = 1.25. This en- 
hanced growth applies in locally overdense regions. Pois- 
son fluctuations from assigning N-p = 1.6 x 10^ particles 
to A*b = 256^ bins generate density fluctuations with a 
standard deviation of ~ ^/N^JN^ « 0.2. The TSC as- 
signment smooths these fluctuations somewhat, but ran- 
dom overdensities 25% or greater still exist in over 1000 
cells (1.6% of the total). Since a region with e « 1.25 is 
already non-linear after two e-foldings (consistent with 
the observed t « 6J7~^), enhanced growth in overdense 
regions plausibly explains the growth of cavities. 

We conflrm this physical explanation for the cavities 
by running flve variations to AB: (1) doubling the spa- 
tial resolution, (2) doubling the number of particles per 
grid cell, (3) seeding a linear mode (an eigenvector) with 
the "cold start" algorithm used for the linear tests in 
YJ, (4) the same linear mode, but the particle density 
distribution is seeded randomly (and thus dominated by 
Poisson fluctuations), and (5) quadratic polynomial in- 
stead of sphne interpolation (see Appendix A of YJ). 
The growth of the root-mean-square of the vertical gas 
velocity for the first three variations is shown in Fig. [6] 
along with the original run BA. Variation (1) [and also 
(4) and (5), not shown] give essentially the same behav- 
ior run BA, which eliminates obvious numerical effects 
(grid resolution and interpolation scheme) as the source 
of cavities.^ Doubling the particle number, variation (2), 
delays the onset of cavitation, as expected with lower am- 
plitude Poisson fiuctuations. Variation (3) suppresses all 
Poission noise, and the "cold" linear mode grows at the 

Variation (5) is interesting because the Poisson density fluctu- 
ations dominate the carefully seeded velocities. 



-2 



-6 



1 1 1 1 1 1 1 1 1 

AB 

■ Double res. 

_ 


1 1 1 1 1 1 1 1 1 1 

Double part. 

Linear mode 












/ 




/ 




/ 




/ 




/ 




/ 




y_ . 


/ 

. . . . 1 . . . . 


^ 

1 .... 1 ... . 



5 10 15 20 

Fig. 6. — Growth of the streaming instability, as measured by the 
amplitude of vertical velocity fluctuations, for different numerical 
approaches to run AB. Only a seeded linear mode (gray line) grows 
with the expected growth rate of s f» 0.1J7 (some curves have 
been wrapped around several times to allow simulations to run 
further). Fast growing cavities occur both at double resolution 
(dotted line) and with twice as many particles per grid cell (dashed 
line). Since cavities are triggered by Poisson density fluctuations, 
explosive growth is delayed with twice as many particles. The 
saturated state is the same for all cases. 

analytic rate, s « O.li?, until non- linear effects finally 
dominate after t = 60J?^^. Perhaps most importantly, 
all approaches lead to the same saturated state, despite 
markedly different routes to turbulence. This speaks to 
the robustness, not just of transient cavitation, but of all 
the non-linear results. 

We also investigated the velocity structure at the onset 
of cavitation. Quadrupolar structures (most prominent 
in the vertical velocity) appear as isolated modes of the 
streaming instability. The length scale of the quadrupo- 
lar distortions did not vary upon doubling the grid reso- 
lution (with a fixed number of particles per grid cell). 

Fortunately, our Poisson noise hypothesis does not pre- 
dict cavitation where it should not occur. Run AC 
(ts — 0.1, e — 3.0) has a fast linear growth rate with 
a relatively weak dependence on the local value of e (see 
Fig.[T]). Accordingly non- linear fluctuations appear uni- 
formly throughout the grid in run AC, instead of cavi- 
tating first in a few spots. The saturated state of run AC 
(shown in the right panel of Fig. [7]) is similar to run AB, 
but with smaller scale fiuctuations. Like run AB, the 
marginally coupled run BB has equal densities of parti- 
cles and gas, but with Tg = 1.0 the rise in growth rates 
across e = 1 was much smoother (see Fig.[T]).^ Since the 
effect of Poisson fiuctuations is weak (an overdensity of 
25% only cuts the growth time by 12% for — 1.0 in- 
stead of halving it for Ts = 0.1) run BB displays orderly 
growth of the dominant linear modes. 

3.3.2. Ts = 0.1, e = 0.2; Weak Overdensities 

Fig. [8] plots the maximum particle density versus time 
for the three tightly coupled simulations. The streaming 

* Fig. 3 of YG confirms this trend, showing that the transition 
is yet sharper for Ts = 0.01 "pebbles." 
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Fig. 7. — The saturated states of run AA (at t = 100017"^) and run AC (at t = 50J7~ ^). Th e range in solids-to-gas ratio is 0.15-0.25 in 
the left plot and 0-15 in the right plot. Run AA (calculated with the two fluid code, see ^3.3.2l for explanation) is dominated by oscillatory 
motion of slightly overdense clumps. The turbulent state of run AC is very much like run AB, but at smaller scales. Also the non-linear 
state of run AC develops simultaneously throughout the grid, unlike the cavitation of run AB. 
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Fig. 8. — Maximum bulk density of solids for the three tightly 
coupled runs (AB, AC and in the insert AA at two different reso- 
lutions). The maximum density is around an order of magnitude 
higher than the average for the e = 1.0,3.0 runs, whereas the tur- 
bulent state of the e = 0.2 run only experiences very mild relative 
overdensities of around 20% at both 256^ and 512^ mesh resolution. 

instability produces particle overdensities of 20 or more 
in runs AB and AC. However the gas-dominated case 
AA has a qualitatively different behavior. Growth sat- 
urates (see left panel of Fig. [7|) in a few growth times, 
igrow = 1/s ~ 42J?~^, as expected. However the par- 
ticle overdensities are very mild, only 20% on average 
(see inset of Fig. [S]) . To test for convergence we ran the 
simulation at both 256^ and 512^ grid points, but the 
qualitative evolution of maximum bulk density of solids 
is unchanged (after a small initial peak in the 512^ run). 
We emphasize that run AA was performed with the 



two-fluid code, not the particle-fluid approach used in 
the other simulations. This choice was necessitated by 
computational cost of long growth times with short wave- 
lengths that restrict the code to small time steps. It is 
tempting to suspect that the weak overdensities in AA 
are a consequence of the two-fluid approach. However 
the hmitation of the pressureless fluid model of sohds 
is that density gradients steepen and shock, causing nu- 
merical instabilities, not that they are stably smoothed. 
To conflrm this we ran two-fluid simulations of case AB 
and obtained the expected (from the particle-fluid run) 
growth of particle density until the code crashed after 
the growth of non-linear overdensities. By contrast AA 
simply never generates large density fluctuations, appar- 
ently since drag feedback on the gas is too weak in the 
small particle, gas-dominated regime. 

3.4. 3-D Simulations 

We have also performed full 3-D simulations of the 
streaming instability with 128'^ grid points and A^p = 
2 X 10^ particles. The linear analysis of YG assumed ax- 
isymmetry, partly for simplicity, but also because modes 
that grow slowly (with s < i7) will be sheared into 
rings in any event. Computer simulations are needed 
to determine whether the saturated state remains az- 
imuthally symmetric in 3-D and how the presence of the 
azimuthal direction affects turbulent properties. Note 
that even axisymmetric linear instabilities can give rise 
to non-axisymmetric parasitic instabilities (e.g. Kclvin- 
Helmholtz instabilities that feed off the channel flow of 
the magnetorotational instability, see I Goodman fc Xul 
IT99I . 

Fig. [9] shows the particle density for runs BA-3D and 
AB-3D in a saturated state. The marginally coupled 
case (BA-3D) maintains a high degree of axisymmetry. 
The radial-vertical plane shows the cascade into sheets 
similar to the 2.5-D case (as seen in Fig. [2]). The quanti- 
tative analysis of turbulent properties (see Tables [2] and 
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Fig. 9. — Saturated streaming turbulence for run BA-3D (e = 0.2, Tb = 1.0, left) and run AB-3D (e = 1.0, Tb = 0.1, right). The boxes are 
oriented with the radial x-axis to the right and slightly up, the azimuthal y-axis to the left and up, and the vertical 2-axis directly up. The 
contours show the particle density at the sides of the simulation box after the streaming turbulence has saturated. The axisymmetry of 
the marginally coupled particles witnesses the smearing effect of Keplerian shear on the relatively long-lived clumps. The tightly coupled 
particles drive rapid fluctuations that develop fully non-axisymmetric density patterns. 



[3]) confirms that BA-3D is very similar to the 2.5-D case. 
The ability to maintain azimuthal symmetry suggests (as 
we will confirm in i i4.3|) that particles reside in clumps 
for longer than an orbital time, so that clumps become 
azimuthally elongated by radial shear. Notice that the 
clump lifetime is not so long that structures appear per- 
fectly axisymmetric. 

The tightly coupled case (AB-3D) on the other hand 
evolves completely non-axisymmetrically. Indeed the 
correlation time of the clumps is short enough that they 
are not significantly elongated by Keplerian shear. Sim- 
ilar to the 2.5-D case, cavities (now fully 3-D and non- 
axisymmetric) developed out of the initial Poisson noise 
in run AB-3D. The saturated state appears to have less 
pronounced clumps than run AB (the fourth panel of 
Fig. [5]). Tables [2] and [3] show that the 3-D turbulence 
indeed has lower velocities (Mach numbers) and weaker 
diffusion, particularly in the vertical direction. It is to 
be expected that turbulent properties in the 3-D runs 
change more for the case that is non-axisymmetric (AB- 
3D) than the case that remains axisymmetric (BA-3D) 
and was already capturing the relevant physics. 

The peak particle densities for the two 3-D runs are 
shown in Fig. 1101 Compared to the density evolution of 
the 2.5-D runs (Figs. Hand [8]) it is evident that BA-3D 
agrees well with BA, whereas AB-3D achieves a some- 
what lower maximum density than AB does. We will 
focus most of our analysis on the 2.5-D runs because 
we could conduct a more systematic study of parameter 
space at higher spatial resolution. The 3-D runs pre- 
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Fig. 10. — Evolution of the maximum bulk density of solids in the 
two 3-D simulations (notice the different time axes). The density of 
AB-3D is somewhat lower than in the 2.5-D case (Fig.|8ll, whereas 
the marginally coupled BA-3D shows good agreement with run BA 
in Fig. |3 

sented here support this choice by giving qualitatively 
(and for BA, fairly quantitatively) similar results to the 
2.5-D runs. 

4. PARTICLE CONCENTRATION 



10 



Johansen & Youdin 



TABLE 3 

Turbulent Transport 



Run 


Ts 


e 










.j-.(turb) 




t:.{NSH) 




AA 


0.1 


0.2 


(1.4 ±6.2) X 10" 


7 


(6.0 ± 262) X 10- 


7 


-2.2 X 10- 


8 


-4.8 X 10- 


7 


AB 


0.1 


1.0 


(4.4 ±0.4) X IQ- 


5 


(2.9 ±0.5) X 10- 


5 


-6.1 X 10- 


5 


-3.1 X 10- 


7 


AC 


0.1 


3.0 


(2.0 ±0.2) X IQ- 


5 


(1.8 ±0.2) X 10- 


5 


-6.0 X 10- 


5 


-5.8 X 10- 


8 


BA 


1.0 


0.2 


(2.2 ±0.6) X 10- 


3 


(1.5 ±0.8) X 10- 


2 


6.7 X 10- 


5 


-1.7 X 10- 


4 


BB 


1.0 


1.0 


(7.6 ±0.7) X 10" 


4 


(1.7 ±0.4) X 10- 


4 


-4.0 X 10- 


5 


-2.0 X 10- 


4 


BC 


1.0 


3.0 


(2.8 ±0.2) X 10" 


4 


(6.2 ±0.9) X 10- 


4 


-1.5 X 10- 


4 


-5.2 X 10- 


5 


AB-3D 


0.1 


1.0 


(1.6 ±0.2) X 10~ 


5 


(2.7 ±0.1) X 10- 


6 


-1.5 X 10- 


5 


-3.1 X 10- 


7 


BA-3D 


1.0 


0.2 


(2.0 ±0.3) X 10- 


3 


(8.2 ±2.5) X 10- 


3 


6.0 X 10- 


5 


-1.7 X 10- 


4 



Note. — Col. (1): Name of run. Col. (2): Friction time. Col. (3): Solids-to-gas 
ratio. Col. (4)-(5): Turbulent diffusion coefficient in units of c^Q~^ (interval indicates one 
standard deviation in each direction). Col. (6): Radial ffux of azimuthal momentum relative 
to NSH state. Col. (7): Radial flux of azimuthal momentum in NSH state. All quantities 
are normalized with standard combinations of 17, Cs and pg. 



The ability of drag forces to concentrate particles via 
the non-linear evolution of the streaming instability is 
now analyzed in detail. This fundamentally important 
process could alter the coUisional evolution of the size 
spectrum of particles, leading to an enhanced growth 
of the average particle radius, or even trigger gravita- 
tional instabilities in the solid component of protoplane- 
tary disks. 

4.1. Gas Does Not Clump 

We emphasize that gas densities remain nearly con- 
stant, despite non- linear particle overdensities in stream- 
ing turbulence. Gas overdensities are < 1% in all runs. 
This validates our use of a constant stopping time, Tf 
(which would otherwise vary with gas density in the Ep- 
stein regime). We note that the linear analysis of YG 
assumed a perfectly incompressible gas. YJ confirmed 
that the linear growth is indeed unaffected by gas com- 
pressibility, which we now see also remains weak in the 
non-linear regime. 

The gas fluctuations are consistent with the small 
Mach numbers in Table [H which are below (but near) 
the scale set by the pressure supported velocity, rjv\^, 
with 7]V]<^/cs — 0.05 in our simulations. Curiously, the 
range in radial Mach numbers is remarkably narrow for 
all the 2.5-D simulations, from 8.7 x 10"^ to 1.2 x 10"^ 
(with the weakly turbulent, two- fluid run AA excluded). 

4.2. Particle Density Distribution 

To get a clear picture of both typical and maximum 
particle overdensities. Fig. [TT] plots the cumulative dis- 
tributions of particle density during the saturated phase 
of the simulations. The distributions measure the frac- 
tion of particles with ambient densities above a given 
value, and are averaged over many snapshots to ensure 
adequate sampling. Particle densities relative to the gas 
are readily obtained by multiplying the x-axis values by 
e = (pp)/(pg). Run BA (r^ = 1.0, e = 0.2) has the 
largest particle overdensities, of nearly 1000, meaning pp 
reaches nearly 200 times the gas density. However since 
run BC (r^ = 1.0, e = 3.0) starts with a particle den- 
sity 15 times larger, it experiences larger peak values of 
Pp/(Pg) ~ 900. Curiously run BB (r, = 1.0, e = 1.0) 
is not an intermediate case but has smaller overdensities 
relative to both particles and gas. 

Particle concentration is more modest during the 
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Fig. 11. — Cumulative particle density distributions. The curves 
show the fraction of particles with an ambient density > pp (the 
dashed line indicates a 10% border between typical and excep- 
tional). Except for run AA [js = 0.1, e = 0.2) the majority of par- 
ticles reside in clumps overdense by a factor of 2—10. A small frac- 
tion of particles experience extreme overdensities of nearly 1000. 

tightly coupled runs. Case AB and AC have very simi- 
lar particle overdensities, with an average (5p w 2-3 and 
a peak 5-p « 30. For case AA (ts — 0.1, e — 0.2) the 
overdensities are negligible. It is remarkable (if a bit 
mysterious) that the e = 0.2 runs give both the strongest 
(BA) and weakest (AA) particle overdensities, depending 
on stopping time! 

Fourier spectra of the particle density are shown in 
Fig. [121 The absolute value of the Fourier amplitudes, 
normalized by the mean bulk density of particles, has 
been averaged over many snapshots during the saturated 
turbulent state of the simulations. Runs BA, BB and 
BC show clear peaks at large scales in agreement with 
the scale of the clumps seen in Figs. [2] and [H whereas the 
power in the tightly coupled runs AB and AC is largely 
isotropic and monotonically decreasing with decreasing 
wave length. Run A A is extremely dominated by the 
very largest scales of the box. 
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Fig. 12. — Power spectra of the bulk particle density, along the x-direction (black line) and the z-direction (gray line). The Fourier 
amplitudes are shown normalized with the mean density of particles in each simulation. Runs BA, BB and BC show clear peaks at large 
scales in agreement with the scale of the clumps seen in Figs. [2] and |3] whereas the power in the tightly coupled runs AB and AC is largely 
isotropic and monotonically decreasing with decreasing wave length. Run AA is extremely top heavy with power almost exclusively at the 
few largest scales of the box (see insert). 

its peak value. Runs BA and BC have the longest 
icon- ~ (6-7) J7^^. Run BA is the best sampled and shows 
a secondary peak past t = 60f2~^ indicating either peri- 
odicity or (more likely) secular changes from the ongoing 
cascade and small clump numbers. Run BB enjoys a 
shorter icon ~ 3.5f2~^, but Cp does not quite drop to 
zero, an indication that a fraction of particles remain in 
dense regions. 

The runs with tighter coupling of — 0.1 (AB and 
AC) had quite short icon- ~ 0.7f2~^ (note the compressed 
time axis in Fig. [T3] for these runs). The short correla- 
tion times are consistent with the less pronounced clump- 
ing and lack of upward cascade when compared to the 
marginally coupled runs. For run AC, Cp is significantly 
negative for t > 20J7^^, indicating that particles avoid 
dense clumps after leaving them. 

It is clear from movies of the simulations that many 
clumps persist longer than icom particularly in the 
marginally coupled Ts = 1.0 runs. The particles that 
make up a clump continuously leak out downstream to 



4.3. Correlation Times 

The residence time of particles in dense clumps affects 
the cosmogonical processes, e.g. gravitational collapse or 
chondrule formation, that might occur therein. For this 
purpose, we measure the time correlation function of the 

ambient density, pp \ experienced by particle i, 



(pW(Op«(t' + 0)-(p«)2, 



(3) 



from snapshots of the particle positions taken every 
At = fl^^ apart. The brackets indicate an average over 
time® and the particles tracked (10% of the total was 
more than sufficient for convergence). Subtraction of the 

mean squared ensures that positive (negative) val- 
ues of Cp{t) correspond to correlation (anticorrelation) , 
respectively. 

Fig. [T3] plots the time correlation function for the sat- 
urated state of the 2.5-D simulations. A characteristic 
correlation time, icorr, is obtained when Cp drops to half 

^ Since averaging is restricted to intervals t apart, the largest 
t considered is never more than half the (non-linearly saturated) 
duration of the simulation. 



We considered defining correlation functions and times only 
for particles initially residing in overdense regions, but equation 
((3)l is a quadratic measure that already favors such regions. The 
simpler, more standard definition is sufficient for our purposes. 
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Fig. 13. — Time correlation functions for particle density indicate 
how long particles reside in dense clumps. Runs BA and BC have 
long correlation times of ~ 7f2~^, while BB has a slightly shorter 
value of 3.51?-^ These Ts = 1.0 runs use black lines and the 
lower black time axis. The tightly coupled runs, AB and AC (gray 
lines and the upper gray compressed time axis), have very short 
correlation times < 117"^. 

the radial drift flow and are replaced with new parti- 
cles drifting in from upstream. The mismatch between 
clump lifetime and density correlation time is evidence 
that the clumps are a dynamical, collective phenomenon 
in the solid component, rather than a persisting density 
enhancement. That situation might change with the in- 
clusion of the self-gravity of the solid particles, as this 
could cause the clumps as a whole to collapse under 
their own weight, fragmenting perhaps into gravitation- 
ally bound objects. We plan to include the self-gravity 
of the particles in a future research project. 

4.4. Energetics of Clumping 

The growth of particle clumps shields solids from the 
full brunt of drag forces, akin to the drafting practiced 
in bicycle pelotons. In YJ §5.1 we show that the rate of 
energy dissipation by drag forces. 



(4) 



is diminished (brought closer to zero) by particle clump- 
ing in the laminar state (at least for tight or marginal 
coupling) . To determine the relevance of this process for 
the saturated turbulent state, Fig.[T4]plots the time evo- 
lution of the energy dissipation rate for the marginally 
coupled runs (black time axis) and two of the tightly cou- 
pled runs (gray time axis). For the tightly coupled runs 
(AB and AC) the dissipation actually becomes stronger 
(more negative) in the saturated state, a consequence of 
increased relative velocities in the turbulent state. Since 
these runs also have significant overdensities, the low- 
ering of l^^dragl is apparently not a necessary condition 
for clumping. The short clump lifetimes (see Fig. fO)) 
is consistent with the inability to reduce dissipation by 
drafting for Tg = 0.1. 

By contrast, all marginally coupled runs (BA, BB 
and BC) show diminished dissipation in the non-linear 
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Fig. 14. — The energy dissipation rate [normalized to 
(pg) {rivj<_)^ n] from drag forces between solids and gas. Marginally 
coupled runs BA, BB and BC (black curves) reduce the dissipa- 
tion rate in the turbulent state, by shielding particles in dense, 
long-lived particle clumps. Tightly coupled runs AB and AC (gray 
curves, which follow the top gray time axis) show increased dis- 
sipation, since particle clumps are too short-lived to allow such 
shielding. 

state, more consistent with the analytic expectations 
from clumps in a laminar flow. The longer correlation 
times and the upward cascade into large clumps exhibit- 
ing bulk motion (particularly in AB, see i)3.2p foster the 
reduction of |5drag|- The resulting particle overdensi- 
ties are significantly larger for these Tg = 1.0 runs (see 
Fig. [11]). Thus diminishing drag dissipation is not re- 
quired to generate particle overdensities, but this draft- 
ing mechanism can augment the growth of dense clumps. 

5. TRANSPORT 

In this section we quantify the effect of the streaming 
turbulence on the radial drift of particles, radial momen- 
tum transport and on the diffusive mixing of solids. 

5.1. Radial Drift 

We initially expected that streaming turbulence would 
reduce the radial migration of particles, due to the pro- 
nounced particle clumping. The laminar drift of particles 
slows as [see YJ equation (7c)] 



,,(NSH) _ 



2r, 



(1 



Pp/Pg) 

2 



for ^ » 1, 



(5) 



(6) 



with increasing particle inertia. The results of the sim- 
ulations are more complicated since turbulent velocity 
ffuctuations produce drift speeds that deviate from local 
equilibrium. 

Table [5] lists average radial drift velocities in the tur- 
bulent state, along with the laminar equilibrium values 
from NSH. Radial drift decreases by about 40% during 
run BA, as is also shown in Fig. [151 For the other Tg — 1.0 
runs, BB displays a modest 15% reduction while BC is 
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Fig. 15. — Evolution of the radial drift speed of solids during 
run BA, averaged over all 1,600,000 particles. Streaming turbu- 
lence slows the influx of solids by 40% below the laminar drift 
speed (dotted line) on average, with signiflcant temporal fluctua- 
tions that correlate with peaks in the maximum bulk density of 
particles (Fig.|4ll. 



unchanged despite significant overdensities. The tightly 
coupled runs show marked increases in drift speeds of 
200% for AB and 300% for AC. Note that BA has the 
fastest laminar drift (due to marginal coupling and low 
solids-to-gas ratio), while AC (followed by AB) have the 
slowest laminar drift (because of tight coupling and large 
particle inertia). While the same ordering of drift speeds 
holds in the turbulent state, the range of speeds for differ- 
ent parameter choices shrinks (i.e. the fastest slow down 
and the slowest speed up). We examine this trend in 
detail below. 

Fig. [TBI shows (with black histograms) the distribution 
of drift velocities, averaged over time in the turbulent 
state, for six different runs. For comparison, the location 
of the equilibrium drift velocity is indicated with a short 
vertical line (labeled NSH). The full gray lines plot the 
average ambient particle density for particles in a given 
velocity bin, and should be compared to the dash-dotted 
gray lines that show the laminar relation between particle 
density and drift velocity (from the inversion of eq. [5]). 
The laminar drift velocities have a finite range from 
for infinite particle densities to the single-particle case, 

wi™'"^ = —r]VK for Ts = 1.0 and wi™'"'' ~ — 0.27/tiK for 
Ts = 0.1. The actual velocity range extends beyond these 
limits in the turbulent state. 

First consider the marginally coupled runs (left column 
of Fig. I16p . The velocity distribution is non-Gaussian 
with a clear negative skewness (velocities drop sharply 
at the right side of the Gaussian, with a more gradual 
decline toward negative velocities). The gray lines show 
the expected trend of slower inward drift for higher ambi- 
ent densities. For particles moving radially outward with 
Vx > 0, the average particle density drops with increasing 
speed. This is reasonable behavior since low density par- 
ticle clumps can more readily be fed angular momentum 



and pushed outward by gas fluctuations.^"'^ The extended 
tails of fast-drifting material at low densities are respon- 
sible for the modest (or non-existent for BC) reduction 
of drift velocities, despite the slowing, or even reversal, 
of motion in overdense regions. 

Now consider the tightly coupled runs in the right col- 
umn of Fig. 1161 The velocity distributions are nearly 
Gaussian and extend well beyond the range of lami- 
nar drift velocities, indicating that turbulent fluctuations 
dominate. The peaks are shifted leftward, which pro- 
duces the higher turbulent drift speeds of Table [2l The 
gray curves plot the astounding fact that overdense re- 
gions drift in faster, a reversal of the laminar trend. This 
is seen in the movie of AB where dense clumps snake their 
way inwards while underdense diffuse material races out 
(the snake patterns are visible in Fig. [5] as well) . 

5.1.1. Effective Drag on Clumps 

The tendency for faster migration of dense clumps for 
Ts = 0.1 can be understood as a consequence of an ef- 
fective, macroscopic drag force acting on the clumps. 
The gas inside the clump is tied to the clump, but 
exterior gas passes freely around the surface, exerting 
an effective drag. If the effective friction time of the 
clump is closer to unity than the original Tg, then the 
dense clump will behave more like a marginally coupled 
solid and drift inward faster. This collective drag effect 
is similar to the plate drag m odel of Ekman layers on 
the surface of particle sub disks (jGoldreich fc Wardlll973t 
iGoodman fc Pindoil[200l . 

We estimate the friction time t^^^'^ of a clump of radius 
^ciump as the time required to encounter its own mass, 
-Wciump, in gas.^^ That gives for 3-D clumps (2-D clumps 
give the same final scaling) 



Jeff) 



clump 



Pp ^clu 



Aw 



(7) 



where pp is the bulk particle density inside the clump 
and Aw is the speed of the clump relative to the gas. 
Multiplying each side by Q yields 



(cff) 



lump ^'^K 



Pg w 



Av 



(8) 



If we dare test this heuristic hypothesis, reading the size 
of the clumpy plateaus from Figs. [5] and [7] and the bulk 
density and speed of the dense clumps relative to the gas 
from Fig.[Tni we find for run AB the values Pp/{pg) « 2.5, 
^ciump/(77'') ~ 0.1, Av w rjVK (the velocity difference 
between high density material and low density material 
in Fig. I16p . corresponding to an effective friction time of 



f2T, 



(cff) 



0.25 Run AC has pp/ (pg) « 6, i?ciump/(?7'') 

.(eff) 



0.05, Av « rjVK, yielding a very similar value of Tj 
0.3. Thus our crude estimates show that the clumps 
couple aerodynamically to the gas more loosely than the 
individual particles do, explaining at least qualitatively 
the faster drift of dense clumps and the increase in drag 
dissipation, two surprising features of the Tg = 0.1 runs. 

In the absence of fluctuations and with an outwardly decreas- 
ing pressure, particles only drift inwards. 

This is the valid criterion for high Reynolds number, turbulent 

drag. 
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Fig. 16. — Histograms of the fraction of particles with a given radial drift velocity Vx in the turbulent state (black curves). The short 
vertical lines (labeled NSH) indicate analytical drift velocities in the initial state with no turbulence or clumping. Marginally coupled (B*) 
runs show a slowing of the net drift speed, whereas tightly coupled runs (A*) produce faster infall. The gray lines (following the right 
y-axes) show the average particle density in each velocity bin. For reference the gray dash-dotted lines plot the laminar particle density vs. 
drift velocity relation. The B* runs display the expected decrease in drift speeds with increasing density, whereas the A* runs (surprisingly) 
follow the opposite trend. 



5.2. Momentum Flux 

The radial flux of orbital momentum, J-c,x = PgUxUy + 
PpWxWy, and its eontribution to disk heating are dis- 
cussed in YJ §5. For laminar flow the drag equilibrium 
between solids and gas gives (equation 18b of YJ) 



(NSH) _ 



•qvK 



.(l + e)2 + ^ 



(9) 



The inward transport of angular momentum follows from 
the the slower rotation of the outgoing gas and the faster 
rotation of the incoming particles. The values for J-c.x 
in the saturated turbulent flow are given in Table [3] 

and are decomposed as J^c,x = ^'c^T^'' + ^c™^''-' 
the laminar value and changes caused by turbulence. If 
the turbulence were driven by orbital shear, which re- 
leases free energy via outward angular momentum trans- 
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port, T ^ would be positive. Instead, most runs have 

•^^c^x^'' < 0, which is physically allowed since work done 
by the global pressure gradient powers streaming turbu- 
lence. Only run BA (and BA-3D) has J^i:'"'"^ > 0, but 
the net angular momentum flux is still inward. Thus in 
all our simulations, angular momentum transport acts 
to take kinetic energy out of the motion, at the rate 
£c — (3/2)i7Jv;.a; < (see YJ §5). As in any shear- 
ing box simulation with (shear) periodic boundary con- 
ditions, momentum fluxes are divergence-free constants, 
which prohibits secular evolution. Global simulations are 
needed to fully investigate the role of "backwards" angu- 
lar momentum transport from streaming turbulence on 
disk evolution. 

5.3. Turbulent Diffusion 

The turbulent mixing of particles is usually modeled 
as a diffusive process in which particle motions are de- 
scribed by a random walk for large length-scales and over 
long time-scales. We provide here best fits to the dif- 
fusion coefficients in the radial and vertical directions. 
Since the motion of particles is very complex, and fur- 
thermore the particles arc not passive contaminants but 
the cause of turbulence, we also test the validity of the 
diffusion approximation. 

We track the deviation of particle positions, Xi{t) and 
Zi(t), from their positions at an initial time to when tur- 
bulence has already developed. Here we are not con- 
cerned with the net radial drift of particles, but the 
spreading of the distribution Xi{t) —Xi(to) (and similarly 
for vertical motions, although no systematic motion over 
long time-scales is expected in this direction). Particles 
are allowed to move greater distances than the box size 
by deconvoluting any particles that were transferred over 
the periodic boundaries by the code. For pure diffusion, 
the distribution tends to a Gaussian with a variance, cr^, 
that grows linearly with time.^'^ An example of how the 
particles spread out with time is shown in Fig. [T7] for 
radial mixing in the BA simulation. The diffusion coef- 
ficients, and Dz, are extracted as 



The best fit diffusion coefficients are listed in Table 
[31 using the standard normalization by c^[2~^. The 
dimensionless diffusion coefficients lie in the interval 
10~^ . . . 10~^, ranging from extremely small up to val- 
ues that are comparable to t he diffusion caused by 
magnetorotational turbulence (jCarballido et al.l l2005t 
[johanscn & Klahr 2005). For the smaller = 0.1 par- 
ticles the diffusion is quite weak, < 5 x 10~^. This is 
because more tightly coupled particles trigger weaker 
turbulence with smaller length scales, a result consis- 
tent with smaller linear growth rates and wavelengths for 
lower Ts (see YG). Run BA (r^ = 1.0, e = 0.2) exhibits 
anomalously large diffusion, especially in the vertical di- 
rection, Dz « 0.01. This is due to the significant bulk 
motion of elongated clumps (discussed in i i3.2l) . 

Since particles are allowed to cross the periodic boundaries, at 
late times different portions of the distribution will overlap. This 
just means that only turbulent scales up to a certain length scale 
are considered, something that should not significantly affect the 
integrity of the measurements. 
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Fig. 17. — Spot the platypus. The histograms plot the radial 
distance that particles in run BA have traveled since the reference 
time of to = 20roj.i,. The curve moves inward due to the radial 
drift, while spreading as a random walk with a Gaussian width a 
that increases as the square root of time. The Gaussians are slightly 
platykurtic, or flat-tailed, due to a population of solid particles that 
experience decreased diffusion in the massive particle clumps seen 
in Fig. [2 

The upper limit for the diffusion appears to be set by 
the characteristic length and velocity scales, r]r and tjvk, 
to he D < rj'^VKr ~ rjc^/n, i.e. D < rj = 5 x 10~^ when 
non-dimcnsionalized. Indeed even the extreme Dz in run 
BA only violates this order of magnitude criterion by a 
factor of three. As a consistency check on the diffusion 
coefficients, D^^z ~ ^w'^ z^corr is obeyed within a factor 
of a few, for random velocities, 5w, from Table [2] (for 
the gas, but particle values are not more than 10% 
different) and the correlation times, ^com from Fig. 1131 

It will be interesting to compare these results to strat- 
ified disk models with self-consistent vertical settling, 
where the relevant parameters are Ts and the solids-to- 
gas surface density ratio (instead of e). 

5.3.1. Validity of the Diffusion Approximation 

We performed several tests of the diffusion approxi- 
mation. The time variation of the diffusion coefficients 
should be small, and especially should lack an overall 
deviation from da^ jdt oc constant. This was true for 
most runs, as indicated by the error bars on the diffusion 
coefficients in Table |3l The two exceptions were again 
run BA, which exhibited large fluctuations due to the in- 
teractions between a few large particle clumps, and the 
two-fluid run AA. This run was seeded with tracer parti- 
cles, following the velocity field of the solid fiuid, in order 
to be able to use the random walk approach to measure 
diffusion. The tracer particles exhibited extremely small 
diffusion with a huge fluctuation interval, an effect of 
the weak non-linear state of run AA where periodic bulk 
motion of a few clumps dominates over random motion 
(see Fig. [71). Particles spread out and gather again in 
a way that is distinctly not a random walk, but over 
longer time-scales the particle distribution still spread 
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out as a Gaussian. The enormous error interval indi- 
cates that the turbulent transport is not like diffusion on 
short time-scales. 

We also tested Gaussianity by measuring the skewness 
and kurtosis of the particle displacement distributions. 
Most runs were fairly Gaussian, except for a modest 
skewness, ~ 10%, in the radial (and not vertical) dis- 
tributions, which is readily explained by the interaction 
of the radial drift flow with clumps. The BA run exhib- 
ited a kurtosis of —0.5, i.e. slightly platykurtic or small- 
tailed (see Fig. [T7l) , consistent with transport influenced 
by bulk motions, and not just a random walk. Model- 
ing turbulent transport as diffusion is under all circum- 
stances only an approximation. Still, the turbulent dif- 
fusion coefficient is a good measure of the time-scale on 
which solid particles are mixed by the turbulent motion. 

6. SUMMARY 

In this paper we have shown that solid particles can 
trigger turbulence in gaseous protoplanetary disks via 
the streaming instability and thereby cause their own 
clumping. We have ignored a number of complicating ef- 
fects. Most critical is perhaps the lack of vertical gravity, 
but we believe it was instructive to see how the stream- 
ing instability evolves in a pure model that has exact 
linear solutions first. We plan to include both vertical 
gravity and the self-gravity of the solids in a future re- 
search project. A distribution of particle sizes and phys- 
ical collisions between particles have also been ignored, 
even though coagulation, fragmentation, and coUisional 
cooling are likely relevant in dense particle clumps. As 
the complex behavior of our simple model system shows, 
significant progress on basic physical processes can be 
made before the "kitchen sink" approach is required. 

The most striking consequence of streaming turbulence 
is the growth of overdense particle clumps without self- 
gravity. This effect was previously seen in the non-linear 
simulations of particle settli n g and K elvin-Helmholtz tur- 
bulence bv I Johansen et al] ( 2006af ). In both that work 
and this one, clumping can be a self-propagating phe- 
nomenon. The increased inertia in dense clumps de- 
creases their drift speeds, creating local "traffic jams." 
We saw this behavior for marginally coupled solids, which 
developed the largest relative overdensities, above 100, 
with an upward cascade to long-lived, vertically elon- 
gated filaments. Marginally coupled solids have the high- 
est radial drift speed and are known to exhibit the most 
pronounced drag-related phenomena, so it is not surpris- 
ing that marginal coupling also gives rise to the most dra- 
matic streaming turbulence. While clumping may not in 
itself explain how to keep large amounts of marginally 
coupled particles at large orbital distances (jWilner et al.l 
[2OOO1) , it does provide a rigorous prediction that the spa- 
tial distribution of such solids will not be smooth, but 
will vary on scales of around one gas scale- height. 

A qualitatively different clumping behavior was seen 
for smaller, more tightly coupled solids. Overdensities 



were lower, in the tens, and clumps were smaller scale 
and short-lived. To extend the analogy, these runs ap- 
peared more like a game of bumper cars than a full-scale 
traffic jam. The biggest surprise was the complete rever- 
sal of the laminar relation between drift speed and par- 
ticle density. Dense clumps actually fell in faster than 
particles in voids for the tightly coupled solids. Our 
heuristic explanation is that robust clumps can with- 
stand turbulent boundary layer ffows that sap their an- 
gular momentum, as in an Ekman layer flow. A similar 
explanation has been applied to the surfaces of particle 
sub-disks as the plate drag an s atz ( Goldrcich & War^ 
[1971 iGoodman fc Pindod 120001 : IWdenschilhna ,2003i ). 
The run with the tightest coupling and lowest initial 
solids-to-gas density ratio, and consequently the lowest 
linear growth rate, developed very meek non-linear den- 
sity fluctuations. Thus non-linear clumping appears to 
require either marginal coupling or a moderately large 
background solids-to-gas ratio of around unity or higher. 
Further studies of the streaming instability for smaller 
particles, such as chondrules, would be interesting, but 
are computationally costly (see 

It is hardly surprising that the solids-to-gas density ra- 
tio strongly affects the non-linear state since the stream- 
ing instability relies on particle feedback to influence 
gas dynamics. However, the sharp transition across 
particle-gas equality is remarkable, especially for tight 
coupling - the weak streaming instabilities in the gas- 
dominated regime becom e explosive onc e the s olids-to- 
gas ratio reaches unity. lYoudin fc Shul (|200^ argued 
that this threshold also sets a limit to the quantity of 
solids that can be stirred by the Kelvin-Helmholz insta- 
bility. Vigorous turbulence (if stronger than particles 
themselves can stir) could prevent the accumulations of 
such high midplane particle densities. There is little 
doubt, however, that dramatic events occur whenever 
particle densities reach that of the gas. 

Planetesimal formation models generally involve either 
high particle densities, in a gravitational collapse sce- 
nario, or efficient coagulation to particle sizes for which 
radial drift is no longer a problem. Both scenarios in- 
volve conditions - high particle densities and/or marginal 
drag force coupling - where streaming instabilities can 
abet further growth towards planetesimals by generating 
overdense clumps in the particle component. 
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